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1  Summary 

What  follows  is  a  report  based  off  of  the  research  done  for  this  grant.  It 
was  done  with  the  purpose  of  utilizing  the  Euler  equations  with  an  artificial 
viscosity  sensor  to  accurately  and  efficiently  model  the  wake  of  an  aircraft 
or  rotor  craft.  This  report  also  utilizes  research  from  two  others  who  worked 
under  this  grant  in  Rachel  Reitz  and  Joshua  Favors  whose  work  consisted 
of  testing  the  resolution  and  boundary  placement  as  well  as  comparative 
analysis  of  two  different  types  of  sensors,  one  based  on  entropy  the  other 
based  on  exploiting  the  modes. 

This  report  focused  on  the  use  of  an  artificial  viscosity  sensor  for  three 
dimensional  application  to  the  wake.  For  the  wake  it  focuses  on  validation 
for  modeling  the  vortex  bursting  element.  Two  changes  were  made  to  the 
sensor  to  both  increase  the  speed  of  the  sensor’s  application  as  well  as  make 
the  sensor  more  robust.  The  sensor  was  made  more  robust  by  allowing  for 
the  sensor  to  sense  on  more  variables  by  changing  the  baseline  decays  ad¬ 
ditive  factor.  In  this  case  the  change  made  allowed  for  sensing  on  velocity 
terms  where  the  velocity  may  not  always  have  a  value  greater  than  zero.  To 
increase  the  speed  the  transition  from  a  nodal  system  to  a  modal  system 
was  augmented  so  that  the  flop  count  could  be  reduced.  These  changes  were 
tested  against  two  standard  problems  and  then  applied  to  the  vortex  bursting 
problem. 
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2  Introduction 

2.1  Motivation 

The  world  is  full  of  highly  difficult  engineering  and  scientific  problems  that 
cannot  be  solved  by  conventional  means  such  as  the  Navier-Stokes  equations. 
This  is  why  experiments  are  developed  to  replicate  the  physics,  analyze  it  and 
unwrap  the  physics  occurring  in  these  problems.  There  is  a  problem  with  a 
pure  experimental  approach  though  in  that  the  cost  surrounding  the  develop¬ 
ment  and  running  of  these  experiments  can  be  very  high.  This  problem  has 
fueled  the  development  of  different  numerical  techniques  to  supplement  these 
experiments  in  solving  the  partial  differential  equations  that  govern  physics. 
This  is  not  to  say  that  any  numerical  simulation  can  replace  an  experiment. 
Their  usefulness  lies  in  their  ability  to  provide  guided  experiments  to  reduce 
the  cost  of  running  multiple  live  experiments.  This  development  has  led  to  a 
major  focus  of  modern  numerical  techniques  to  solve  the  partial  differential 
equations  that  govern  the  physics  of  these  problems. 

The  desire  to  always  further  the  study  of  numerical  techniques  has  led 
to  many  vast  improvements  in  both  their  implementation  and  the  types  of 
solvers  available.  A  major  development  was  moving  from  the  single  point 
finite  difference  (FD)  method  into  the  finite  volume  (FV)  and  finite  elements 
(FE)  methods  which  utilize  elements.  The  goal  of  all  these  methods  was 
centered  on  being  accurate  to  the  physics  of  the  problem  while  also  making 
itself  faster  than  previous  iterations  allowed.  To  further  improve  upon  these 
methods  accuracy  a  4th  method  was  developed,  the  discontinuous  Galerkin 
(DG)  method  .  This  method  has  been  introduced  relatively  recently  com¬ 
pared  to  the  other  methods.  This  causes  a  necessity  to  improve  the  scheme 
and  make  it  more  robust  and  efficient  in  order  to  remove  any  qualms  about 
the  method.  This  is  very  common  with  numerical  techniques  whenever  they 
are  developed. 

This  report  utilizes  this  DG  method  to  model  features  of  wake  flow  and 
improves  upon  an  existing  sensor  to  make  it  more  robust  for  a  multitude  of 
problems.  The  eventual  goal  is  to  eventually  develop  a  fully  realized  wake 
modeling  code  utilizing  these  techniques.  The  simulation  of  a  full  wake  is 
highly  important  for  all  types  of  engineering  problems  spanning  from  airport 
takeoffs  and  landings  to  the  development  of  car  drafting  techniques.  There 
is  a  reason  these  problems  are  so  crucial  and  interesting  though  and  that  is 
the  problems  they  present  in  order  to  be  correctly  modeled.  A  very  large 
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problem  is  wrapped  in  the  scale  of  the  problem.  It  is  very  difficult  to  try  and 
capture  all  of  the  physics  of  the  problem  through  a  correctly  sized  domain 
because  it  will  always  lead  to  a  preference  being  placed  on  either  the  small 
or  large  scale  features.  This  is  why  not  only  the  DG  method  can  be  a  vast 
improvement  but  the  adaptive  DG  allows  for  an  even  better  solution.  For 
the  scope  of  this  report  certain  key  features  will  be  examined. 

A  common  feature  in  the  evolution  of  a  wake  is  the  sudden  change  in  a 
vortex  core  with  no  change  in  its  circulation  due  to  a  vortex  bursting.  This  is 
a  well  established  phenomena  and  the  problem  has  been  undertaken  by  many 
modern  day  codes  utilizing  both  the  large  eddy  simulation  (LES)  and  direct 
numerical  simulation  (DNS)  .  These  techniques  are  very  useful  but  also  are 
very  intricate  in  how  they  work  and  thus  require  a  large  amount  of  time  to 
be  run.  By  studying  and  comparing  the  data  developed  from  this  current 
configuration  to  the  ones  already  on  hie  the  goal  is  established  to  prove  this 
method’s  effectiveness. 

2.2  Background 

Since  this  work  is  based  on  a  key  aspect  of  wake  evolution  a  brief  history  of 
the  development  of  wake  analysis  theory  and  its  applications  follows.  This  is 
then  followed  by  a  general  overview  of  the  history  of  the  DG  method. 

2.2.1  A  Brief  History  of  Wake  Theory 

The  concept  of  aircraft  wakes  has  been  an  important  aspect  of  aviation  for  a 
long  time.  The  wake  has  strong  impacts  on  a  multitude  of  situations  such  as 
aircraft  taxiing  and  deployment  of  rotocraft  on  aircraft  carriers.  The  decay 
of  vortices  throughout  the  wake  plays  one  of  the  most  important  roles  in 
this  research  especially  at  airports.  On  a  grander  scale  there  is  a  hope  that 
with  the  development  of  a  better  understanding  of  vortices,  the  same  will 
be  said  of  turbulence.  Turbulence  remains  one  of  the  greatest  mysteries  of 
the  modern  era  while  also  being  one  of  the  most  divisive  topics  to  discuss. 
The  exact  definition  of  turbulence  is  not  something  that  is  even  unilaterally 
agreed  with. 

The  analysis  of  aircraft  wakes  has  been  undertaken  for  as  long  as  aircraft 
have  existed.  Scientists  such  as  Crow  [6]  spent  a  lot  of  time  developing 
analytical  theories  and  models  in  the  1970’s.  Crow’s  work  was  about  an 
understanding  of  the  vortices  contained  in  the  wake  of  an  aircraft  being 
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dispersed  through  natural  means.  Crow  in  fact  developed  the  theory  behind 
a  common  instability  that  has  been  well  validated  by  todays  experiments  and 
simulations.  This  instability  today  is  proclaimed  as  the  Crow  instability.  In 
the  time  between  1970  and  the  mid  1990’s  wake  modeling  had  advanced 
from  a  dependence  on  different  analytical  models  to  accepting  the  use  of 
computational  fluid  dynamics  models  to  supplement  some  of  the  research. 
The  scientists  Lewellen  and  Lewellen  [13]  in  1996  developed  one  of  the  first 
extensive  three  dimensional  models  of  vortex  wakes  using  a  custom  version 
of  LES  where  they  focused  on  different  aspects  of  the  gas  dynamics  and 
chemistry.  They  specifically  focused  on  the  chemistry  of  the  engine  wake 
from  the  aircraft  and  it’s  interaction  with  the  atmosphere. 

The  work  of  Crow  and  others  was  then  brought  together  by  Spalart  [21] 
in  1998  who  provided  an  encompassing  survey  of  the  knowledge  about  wakes 
to  that  point.  He  discusses  a  multitude  of  issues  concerning  the  propaga¬ 
tion  of  wakes  from  aircraft  and  different  aspects  of  analyzing  them  through 
experimental  setups  at  airports.  He  makes  special  note  of  the  struggles  in 
developing  a  model  due  to  the  limitations  of  the  assumptions  needed  in  CFD 
of  that  era.  He  also  helped  develop  some  of  the  base  theory  concerning  the 
vortex  bursting,  challenging  a  preconceived  notion  that  bursting  was  the  de¬ 
struction  of  vortices,  which  holds  no  weight  when  examining  the  argument 
from  the  standpoint  of  circulation. 

These  issues  are  still  being  examined  and  resolved  in  today’s  climate  es¬ 
pecially  concerning  CFD  and  experimental.  Multiple  programs  have  been 
established  such  as  the  Aircraft  Vortex  Spacing  System  (AVOSS)  used  by 
NASA  to  study  these  vortices  and  establish  more  scientifically  based  param¬ 
eters  especially  with  the  influence  of  the  atmosphere.  Gerz  et.  al.  [8]  utilizes 
the  data  from  AVOSS  along  with  other  similar  agencies  discusses  different 
ways  to  discuss  further  change  that  can  be  applied  to  airports.  Their  goal 
was  to  propose  ways  to  increase  the  efficiency  of  airports  by  utilizing  known 
data  about  wakes  in  order  to  increase  the  number  of  commercial  aircraft  that 
can  take  off  and  land  over  the  course  of  a  day. 

A  common  element  of  all  the  research  is  the  breakdown  of  vortices  which 
is  explicitly  examined  by  Holzapfelfll]  et.  al.  This  group  examined  the  ef¬ 
fects  of  primary  and  secondary  vortices,  which  are  created  through  moderate 
turbulence  from  the  ambient  flow,  have  on  each  other  in  the  case  of  decay 
through  numerical  simulations.  It  is  also  examined  by  Moet  et.  al.  [16] .  Their 
research  is  used  as  the  basis  for  the  vortex  bursting  validation  undertaken  in 
this  research. 
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2.3  Development  of  Numerical  Methods 

In  engineering,  and  science  in  general,  the  development  of  suitable  numerical 
techniques  has  been  a  staple  of  research  for  many  years.  Common  methods 
have  been  used  with  the  FD,  FV,  and  FE  methods.  In  terms  of  modern 
CFD  the  majority  of  codes  use  second  order  accurate  FV  [23]  methods  or 
FD  methods.  While  these  methods  are  useful  they  do  not  fulfill  the  long  term 
goals  of  the  numerical  sciences.  They  are  difficult  to  increase  the  order  of 
accuracy  as  the  mesh  complexity  is  increased  moving  from  pure  Cartiesian 
to  non-uniform  meshes  with  triangular  elements.  They  are  also  expensive 
computationally  for  the  spatial  discretization  of  problems.  This  is  not  a 
problem  with  the  Finite  Element  method  but  it  has  major  flaws  in  that  for 
each  time  step  the  entire  matrix  is  solved  at  once  creating  inconsistencies 
when  applied  to  the  physics  of  a  fluid  mechanics  problem  where  information 
changes  with  the  flow. 

These  desires  led  to  the  development  of  a  new  method  of  solving  PDE’s 
that  gave  the  flexibility  required  for  solving  two  dimensional  and  three  di¬ 
mensional  problems  while  also  allowing  for  a  higher  order  of  accuracy  akin 
to  FE  codes.  This  new  method  was  called  the  DG  method.  This  method 
was  first  proposed  as  a  new  way  to  solve  the  neutron  transport  equation 
by  Reed,  etc.  [18].  The  method  slowly  grew  in  its  different  applications 
to  different  fields  but  it  was  not  till  the  mid  1990’s  that  this  method  ex¬ 
ploded  into  applications  to  CFD.  Advances  were  made  in  the  development 
of  Runge-Kutta  DG  methods  that  provided  an  explicit  and  cheap  method  to 
use  in  the  solution  of  problems.  In  1996  [2]  and  1997  [3]  Bassi  and  Rebay 
published  two  papers  creating  the  local  DG  (LDG)  method  and  applying  it 
to  the  Euler  and  Navier-Stokes  equations  respectively.  This  method  further 
improved  the  previous  results  because  in  the  previous  DG  methods  only  the 
time  discretization  was  discontinuous  but  now  both  space  and  time  were. 

The  LDG  method  quickly  became  one  of  the  most  useful  adaptations 
of  the  method.  Cockburn  and  Shu  [4]  analyzed  this  method  and  showed 
the  great  promise  that  comes  from  its  inherent  flexibility.  It’s  discontinuous 
nature  throughout  allows  for  easy  adaptations  to  parallel  processing  which 
would  allow  for  the  code  to  become  much  faster.  They  would  later  compare 
this  method  against  versions  of  the  DG  method[5].  Here  they  found  the  LDG 
method  had  inherent  advantages  in  adapting  to  not  only  elliptical  PDE’s 
but  also  hyperbolic  and  parabolic  PDE’s  allowing  for  a  more  general  use 
of  this  method  to  problems  even  outside  of  the  realm  of  CFD.  This  was 
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further  explored  by  Yan  and  Osher  [24]  in  the  solving  of  the  Hamilton- Jacobi 
problems.  This  was  then  taken  further  by  Kevin  Albarado  when  he  purposed 
the  DG  method  to  solve  the  Hamiltonian  for  the  burning  of  a  two  dimensional 
start  solid  rocket  grainfl]. 

2.4  Objective 

The  purpose  of  this  report  is  to  further  the  development  of  a  robust  DG 
scheme  based  on  the  Euler  equations  and  implementing  an  artificial  viscos¬ 
ity  sensor  that  has  been  augmented  to  make  it  more  applicable  to  different 
this  problem.  Artificial  viscosity  is  a  necessity  in  that  even  though  the  Eu¬ 
ler  equations  are  inherently  invisicid,  the  numerical  equations  solved  have  a 
viscous  term  when  backed  into  their  conservative  forms.  This  is  normally 
taken  care  of  in  the  inherent  dissipation  in  the  numerical  scheme  but  with 
the  high  order  of  the  DG  method  this  needs  to  be  readdressed.  This  artifi¬ 
cial  viscosity  sensor  will  be  altered  and  improved  upon  in  order  to  allow  for 
more  versatility  and  robustness  on  a  range  of  problems  from  shocks  to  vortex 
interactions  while  also  increasing  the  speed  which  is  always  important  when 
working  with  numerical  simulations  and  especially  when  these  problems  take 
place  in  three  dimensional  space.  Different  test  problems  will  be  examined  in 
order  to  accurately  show  the  effects  of  these  two  changes  to  these  important 
factors. 

After  the  sensor  has  been  tuned  for  optimal  use  on  three  dimensional 
problems  the  vortex  bursting  problem  will  be  solved.  This  will  be  done  by 
using  a  complex  method  for  developing  the  initial  condition  that  utilizes  a 
pseudo-steady  state  problem.  The  mesh  size  and  polynomial  order  will  be 
varied  to  show  the  effects  on  the  full  problem.  The  helical  instability  will 
be  tripped  as  well  to  show  a  more  physically  accurate  problem  as  well  as  to 
show  that  the  sensor  can  account  for  artificial  viscosity  where  needed  but  not 
overload  the  problem  and  ’’wash  out”  useful  and  accurate  physics  from  the 
problem. 


3  Discontinuous  Galerkin  Method 

What  follows  is  the  development  of  the  DG  method  and  a  process  to  apply  it 
to  a  simple  one  dimensional  problem.  This  analysis  follows  a  procedure  done 
by  Hesthaven  et.  al.  [10].  The  DG  method  is  a  hybrid  of  sorts  of  the  FV  and 
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FE  methods  in  order  to  gain  the  distinct  advantages  of  each  method.  This 
allows  the  new  method  to  have  the  geometric  flexibility  of  a  FV  code  while 
allowing  for  much  more  flexible  method  of  creating  a  higher  order  of  accuracy 
spatially.  This  is  a  major  problem  in  many  current  CFD  codes  running  FV 
solvers.  The  simple  one  dimensional  transport  equation  will  be  used  as  a 
tool  for  the  application  of  this  method  where  u  represents  the  velocity  of  the 
wave. 


du 

dt 


9f(u) 

dx 


=  0 


(1) 


The  first  part  of  any  numerical  scheme  is  to  take  the  domain  of  the  prob¬ 
lem  (Q)  and  discretize  it  into  K  number  of  cells  for  a  new  domain  (f\)  . 
The  point  of  discretization  is  to  break  a  domain  into  something  that  can  be 
solved  computationally.  Each  of  these  cells  attaches  to  surrounding  cells  at 
an  interface  that  overlaps.  This  is  an  important  feature  for  the  creation  of 
a  flux  that  occurs  between  the  two  cells  which  will  be  discussed  later  in  this 
section.  An  example  of  how  the  discretized  domain  is  represented  is  shown 
in  Figure  1. 


D 


k-1 


D 


D 


k+1 


-J< 

k+1 

x. 


Figure  1:  ID  discretization  of  a  domain  into  three  elements 


The  global  true  solution  u  is  then  first  equated  to  the  computational 
solution  and  then  broken  up  into  each  cell  such  that 

K 


u(x,  t )  ~  uh(x,  t)  =  0  ukh(x,  t )  (2) 

k= 1 


where  there  is  now  a  state  for  each  cell  that  combines  into  an  approximation 
with  the  global  state.  A  high  order  local  basis  function  is  then  applied  to 
each  cell’s  state  in  order  to  create  a  nodal  state 
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(3) 


N 

40M)  =  ^Uj(t)bfe) 

3=0 

The  residual  7 Zh  is  then  created  and  with  this  method  it  has  to  vanish  for 
each  cell  in  a  Galerkin  sense  where  the  basis  function  forms  the  weighting 
function.  This  is  how  this  method  differs  from  a  traditional  Galerkin  FE 
method  where  the  residual  vanishes  globally.  This  is  also  what  allows  the 
method  to  be  more  open  to  a  parallelization  which  will  be  discussed  later. 

=  M  +  M  (4) 

[  b)nh dx=0  j  =  0, . ,N  (5) 

J  Dk 

Here  N  is  the  order  of  the  polynomial  and  for  the  one  dimensional  test  case 
the  number  of  points  per  cell  is  equal  to  one  more  than  the  order  (Np  —  N+l) 

With  the  distinction  that  each  cell  is  unique  the  boundaries  of  each  cell 
can  have  multiple  solutions  depending  on  which  cell  is  being  examined  which 
leads  to  Figure  2. 


D 


k-1 


D 


k+1 


D 


Figure  2:  Cells  having  different  boundary  values 


They  can  then  communicate  between  cells  through  a  numerical  flux  func¬ 
tion  which  in  terms  of  CFD  is  physics  based  with  examples  of  Lax- Friedrichs 
or  AUSM+. 

With  the  cells  decoupled  the  original  function  was  then  solved  by  an 
integration  by  parts  using  the  divergence  theorem  to  obtain 


dukh 

dt 


dx  + 


bkdfh 


dx 


dx 


0 
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(6) 


where  n\ i  is  the  element  normal.  Now  even  though  the  cells  are  independent 
of  each  other  there  is  still  a  piecewise  continuous  element  of  the  solution 
concerning  the  flux.  As  previously  mentioned  this  is  a  problem  because  all 
the  cells  needs  to  be  discontinuous.  By  using  the  previously  mentioned  flux 
function  the  previous  equation  now  allows  for  the  discontinuity  at  the  cell 
edge. 

xedD1':  -»/*(«»  (9) 

The  flux  is  designed  such  that  the  flux  coming  from  one  cell  is  equal  to  the 
negative  of  the  flux  leaving  the  adjacent  cell  at  the  same  boundary. 


This  gives  Np  *  K  unknowns  and  equations  which  can  be  solved  cell  by 
cell.  Here  lies  the  main  question  of  the  DG  method  which  is  to  find  u\  G  Pat 
such  that  bi  G  Pat.  Where  P^v  stands  for  a  function  space  on  the  standard 
interval  of  order  less  than  or  equal  to  N  and  vanishes  everywhere  else.  This 
leads  to  a  need  for  a  strong  basis  function. 

3.1  Polynomials  for  Cells 

3.1.1  Cell  Mapping 

An  important  aspect  of  a  Galerkin  model  is  a  strong  basis  function  to  work 
with.  Polynomials  are  widely  used  here  because  of  their  relative  simplicity 
over  Fourier  modes.  There  is  an  exception  in  the  case  of  spectral  methods 
due  to  Fourier  models  providing  a  simpler  model  to  work  with.  They  also 
follow  closely  with  how  people  tend  to  think  about  data  in  the  first  place.  A 
standard  interval,  x  — *  £  G  (—1,1),  is  used  to  develop  a  mapping  function 
which  can  be  extrapolated  to  any  standard  cartesian  cell. 

x  G  Dk  :  x  =  Mk({)  (12) 
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(13) 


Mk( {)  =  x\  +  (XJ  _  4) 

3.2  Basis  Functions 

The  basis  functions  were  correlated  with  the  state  such  that 

N  N 

<4«)  =  £  fiiMQ  =  £  (14) 

m= 0  n=0 

Where  uk  corresponds  to  the  modal  state  and  uk  the  nodal  state.  This 
is  particularly  useful  because  both  nodal  and  modal  analysis  have  different 
strengths  for  analysis.  This  ability  to  move  between  the  two  gives  this  method 
another  advantage  over  traditional  FV  methods. 

3.2.1  Modal  Analysis 

A  traditional  thought  in  the  creation  of  a  modal  analysis  is  the  development 
of  a  monomial  based  on  the  given  polynomial  order  to  set  the  basis  function. 
Then  using  an  L2-projection  to  recover  the  state  uk%  leads  to 


Mu  =  u 


(15) 


The  problem  is  that  this  system  fails  clue  to  the  matrix  M  being  very  poorly 
conditioned  and  thus  the  modal  state  cannot  be  recovered  from  an  inverse  of 
M .  This  led  to  the  idea  of  using  a  L2  based  Gram-Schmidt  orthogonalization 
in  order  to  develop  a  basis  that  was  orthonormal. 

\/ Pm+l'Rm+l  (0  (£  (^m)’^’m(0  \J , 1  (0  ,  . 

!—  (16) 
7T_i  =0,  7T0  =  1/VA) 


This  resulted  in  the  choice  of  the  Legendre  polynomials  which  is  found  by 
selecting  certain  values  for  am  and  j3m.  The  Legendre  coefficients  are  listed 
as 


f  2  if  m  =  0 

[  1/(4  —  m~2)  ifm  >  1 


(17) 


An  important  property  of  these  legendre  polynomials  is  that 


^(0^(0^  =  k 


1  if  i  =  j 
0  if  i  ±  j 
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(18) 


A  modal  basis  is  not  as  useful  though  in  a  direct  analysis  of  the  data  due 
to  the  extra  calculations  that  will  need  to  be  undertaken  to  translate  the 
boundary  conditions  and  initial  conditions  to  a  modal  system.  It  is  though 
a  great  system  that  works  well  for  the  development  of  the  sensor  which  will 
be  detailed  later.  A  transformation  matrix  is  used  to  move  between  nodal 
and  modal  values  for  the  state. 

3.2.2  Nodal  Analysis 

An  ideal  nodal  basis  function  is  the  Lagrange  equation  due  to  it  guarantee¬ 
ing  nodal  accuracy  for  any  quadrature  which  could  be  used.  The  Lagrange 
polynomial  is  created  through 


N  P-P 

4«)=  n  aa  (i9> 

•  n  •  /  S Si 

and  contains  the  following  powerful  property 

4(a=%={  j  (20) 

With  a  powerful  nodal  representation  the  quadrature  needs  to  be  es¬ 
tablished  based  around  the  Legendre  polynomials  needs.  This  leads  to  the 
Legendre- Gauss-Lobatto  quadrature  detailed  below. 

3.2.3  Numerical  Quadrature 

Quadrature  is  an  important  decision  in  the  development  of  any  numerical 
technique.  This  especially  holds  true  if  one  wants  to  perform  numerical  inte¬ 
gration  which  is  based  on  the  distribution  of  the  quadrature  and  is  performed 
as  follows. 

/  <21) 

J~l  j= 1 

Here  c Oj  represents  the  weights  associated  with  a  particular  quadrature  Since 
a  nodal  basis  has  been  already  chosen  for  the  development  of  this  method 
Gaussian  quadrature  seems  the  most  logical  due  it  being  an  exact  integration 
technique  for  polynomials. 
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In  the  development  of  this  method  the  Gauss-Lobatto  rule  was  applied 
since  the  Gauss  rule  does  not  include  the  boundary  points  thus  a  Jacobi- 
Lobato  matrix  was  created. 


1L  — 
^r.J-9  - 


P+2  /73* 

LV+p+mp+i 


Jp+ 1  \/i^p+ieP+i 


a 


p+i 


where  Jp  represents  the  Jacobi  matrix 


«o  \fWi 
\fW\  Oil 

0  . 


0 

0 

0 


CXp—l 


(22) 


(23) 


The  orthonormal  property  of  the  modal  analysis  needs  to  be  kept  in  place  so 
the  application  of  the  legendre  values  for  a  and  f3  were  used.  This  leads  to  the 
weights  and  node  locations  being  found  from  the  eigenvalues  and  eigenvectors 
for  the  Jacobi-Lobatto  matrix  where 


&  =  Af,  Ui  =  P o(^p))2,  *  =  0,  ...,p+l  (24) 

here  A  represents  the  eigenvalues  and  v  represent  the  eigenvectors. 

The  application  of  the  Lagrange  basis  to  these  points  that  satisfy  the 
Legendre  polynomial  is  not  a  problem  and  thus  the  two  methods  now  have 
a  standard  quadrature. 


3.2.4  Combining  the  Two  Analysis 

The  transformation  matrix  between  the  modal  and  nodal  states  was  some¬ 
thing  needed  and  was  developed  through  the  application  of  the  Lagrange 
polynomials  to  the  same  points  that  the  Legnedre  polynomials  were  built 
around.  An  example  of  how  these  two  polynomials  look  up  to  a  5th  order 
can  be  seen  in  Figure  3. 
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(a)  Modal  basis  with  7rm(£)  £  Pat 


(b)  Nodal  basis  with  £„(£)  £  Pat 


Figure  3:  Basis  Function  of  P 


N 


By  applying  the  legendre  polynomial  as  the  basis  function  for  u  it  is  found 


that 


Nv 


(25) 


771=1 


Where  Pm-i  represents  the  legendre  polynomials  which  leads  to 


Vu  =  u  (26) 

where  V  represents  a  generalized  Vandermonde  matrix.  Now  both  the  Van¬ 
dermonde  matrix  and  its  derivatives  can  be  related  to  the  legendre  polyno¬ 
mials  and  their  derivatives  such  that 

U,  '  (Vc),,  '  tt'(&)  (27) 

This  helps  to  create  a  bridge  between  Lagrange  and  Legendre  polynomials 
so  that  the  state  can  be  altered  between  a  nodal  or  modal  value.  This  allows 
the  DG  method  to  take  advantage  of  the  unique  properties  provided  by  modal 
and  nodal  analysis  in  the  orthonormality  of  the  modal  function  with  the 
simplicity  of  the  nodal  system. 

N  N  N 

vK)  =  Eb(«)  =  IWVar)  •*=>  m  =  (28) 

i=0  i= 0  j= 0 
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N  N  N  N 

^/i(0  ^  ^  ^  '  Ujn^n  (Q  ''  ''  ^  ^  ^imMm  ^  ^  ArAn  (29) 

m=0  n=0  m=0  n=0 

The  Legendre-Gauss-Lobatto  quadrature  leads  to  a  well-conditioned  ma¬ 
trix  in  the  Vandermonde  matrix.  With  the  basis  established  the  next  step 
entails  the  manipulation  of  Equation  11  to  be  of  an  operator  form  that  can 
be  solved  numerically. 

3.3  Operator  Form 

Returning  now  to  the  initial  one  dimensional  problem  Equation  1  and  its 
discretized  form 

L b‘°4dx - L + L b‘r (“‘"' +; ”*■) ds = °  (30) 

In  order  to  solve  the  problem  computationally  the  discretized  form  of  the 
equations  needs  to  be  broken  down  into  a  matrix  form  based  on  the  applica¬ 
tion  of  the  basis  function  to  the  state.  The  flux  is  then  broken  up  the  same 
way  the  state  is  based  on  the  basis  function 

N 

=  (3i) 

1=0 

The  time  derivative  term  is  broken  down  through  a  coordinate  transfor¬ 
mation  between  x  and  £.  This  transformation  is  based  on 

x  G  Dk  :  dx  =  Jk df ,  Jk  =  X»~X"  (32) 

and  is  then  inserted  into  the  time  derivative  term  such  that 
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The  same  method  is  then  applied  to  the  spatial  derivative  term  in  order 
to  break  it  down  into  a  derivative  matrix  which  can  be  applied  ot  the  flux 
term 


dx 


N 

£ 

3=0 

N 

£ 

3=0 


dbKxJUk 


I  Dk  dx 

f1  MO 

l-i  ^ 


b’Ax)dx  f: 


bj(Z)<%)f} 


(34) 


N 


=  Z(s%ff 


The  boundary  term  is  trickier  in  that  specific  nodes  based  on  faces  and 
edges  have  to  be  used  in  the  application  of  the  flux.  This  does  not  encompass 
the  global  flux  running  through  the  system. 


dDk 


bif*(ukh  ,v%+;n% 


l/c  s*  (  k. —  k,-\~  k, — 

bif  (  Uh  ,uh  Vnx 

N 


+ 


7  k  (  k, —  k,-\-  k, — 
bif  ( uh  >uh  ;  nx 


N 


=  if(4)  £  +  ft f(4)  £  (K)X(4) 

m= 0  n= 0 

N  _ _  N  _ 

=  v  (6‘  (4)  bt  P3)  (f;)m  +  ]T  (6?  (x{)  (4))  (/; 

m= 0  n=0 

N  _ k  N  _ 

=  Y,  (fc(-i)U-i))  (f:)m  +  £  (M1)M1))  (A- 

m=0  n— 0 

*  _ fc  *  _ k 

=  ( fa)m  +  X!  (fb)n 


k  (35) 


k 

I  b)  n 


m= 0 


71=0 


This  leads  to  the  original  discretized  form  of  the  equation  going  from 


duf; 


- 1  ^  + 1  vr  ( «;r ,  «;r ;  =  o  (36) 


<%? 


at 


I  Dk 


dx 


dDk 
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to 


TV 


TV 


TV 


TV 


j= 0 


dt 


E  (W,  + E  (iJ.™  (/.-)„ + E  («)„  <37> 


»n=0 


m=0 


11=0 


Where  there  are  now  mass,  stiffness,  and  lifting  matrices  for  the  development 
of  the  method.  This  can  further  be  reduced  into  a  matrix  vector  form. 


r]  - - -  —  h  - - -  —  I 

JkM--STfk  +  La(f*)  +La{fi ) 


(38) 


3.4  Elementwise  Operator 

As  previously  mentioned  the  nodal  method  has  been  chosen  in  this  research 
due  it  having  an  advantage  in  dealing  with  boundary  data  as  well  as  non¬ 
linear  flux  functions.  But,  there  is  no  reason  not  to  exploit  some  of  the 
greater  properties  of  the  modal  transformation  in  its  orthonormal  basis  and 
its  hierarchal  construction. 

The  ability  of  the  Vandermonde  matrix  to  manipulate  the  state  between 
the  two  basis  is  a  huge  advantage  in  the  development  of  the  mass,  stiffness, 
and  lifting  matrices.  This  can  be  seen  in  the  development  of  the  mass  matrix 

Mi,  =  /  '  m^M) 

1  TV  TV 

=  /  EET)i.'»(flE(r\'.(a« 

^  1  m= 0  n=0 

TV  TV  /  ,1  \ 

=  EET"T).m  /  (y-% 

m=0  n= 0  W  ^  ' 

TV  TV 

EEVly-fr'), 

m= 0  n=0 
N  N 

EEVyr1),  <=>  M  =  (vvT)-' 

m= 0 n= 0 

Then  with  the  fore  knowledge  that  there  are  recursion-based  formulas 
for  Legendre  derivatives,  the  creation  of  a  derivative  matrix  can  easily  be 
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accomplished  from  the  derivative  of  the  Jacobi  matrices. 


fe) 

N 

m= 0 
N 

=  E  ( v~T)im  (u)™ 

m=0 

N 

=  E  ^  B  = 

m=0 


(40) 


The  stiffness  matrix  for  a  nodal  basis  is  built  by 

sa  =  J  me'jim  ^  =  /W«*  (41) 

which  when  applied  to  the  mass  matrix  and  derivative  matrix  gives 


N 

Y  MinDnj 

n= 0 


E  (//fK)4K)d«)  e'j  «„) 

<>tt)  (E©4tt))d« 


n=0 

/: 

/■1 


mi 

/■l 


0 


,7-1 


5  =  MD 


(42) 


There  is  a  powerful  property  here  due  to  the  need  of  the  transpose  of 
the  stiffness  matrix.  By  applying  the  transpose  and  through  manipulation 
using  the  Vandermonde  matrix  the  development  of  the  stiffness  matrix  can 
be  found  through  just  the  Vandermonde  matrix  and  tis  derivatives  in  the 
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necessary  directions. 

S  =  MD  <=>  ST  =  DtMt 

=  (v^)t  (v-Tv-r 
=  (v^v-Y  (v^v-1)  1  j 

=  {Vt:V-l)T  (VV^-1 

The  development  of  the  lifting  Matrices  are  difficult  though  because  they 
have  to  extrapolate  the  flux  from  the  edge  of  a  cell  to  the  rest  of  a  cell.  The 
suffix  a  and  b  will  be  applied  ot  the  Lifting  matrices  to  denote  which  the  cell 
walls  in  any  direction. 


(£«)«  =  «.(-%-(!) 
=  U  te>)  «,  «<>) 

—  OQjOoj 


M  =  mm 

=  Pi  (€n)  Pj  (€n) 

=  OjyiSisrj 


(44) 


The  numerical  flux  is  then  defined  by  the  same  procedure  to  each  cell  wall 


(/:)‘ =E  (/:)/,■(-!) 

*  ___t 

(A*)‘  =  E  ( WiW 

3=0 

3=0 

N  _ k 

=  52  (KWh') 

N _ k 

=  E(AViK» 

3=0 

3=0 

=  E  (®X 

=  E(T)y« 

3=0 

3=0 

=  (fa) 0 

=  Ub)N 

(45) 


This  all  leads  to  the  combination  of  the  lifting  matrices  and  cell  fluxes 
such  that. 
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(48) 


v  (Lb)ij  u;)j  =  E  w^wij  =  V  ma% 

2=0  2=0  2=0 


These  new  matrices  are  then  plugged  into  the  local  matrix  form  in  Equa¬ 
tion  38  giving  a  complete  system  to  create  the  spatial  derivatives  for  the 
problem. 

4  Implementation  of  Euler  Equations 

For  this  project  the  Euler  equations  were  chosen  due  to  the  large  scale  eddies 
being  the  important  factors.  Therefore  the  small  scale  turbulent  structures 
and  the  turbulence  of  the  problem  is  not  as  important  especially  since  any 
structure  that  needed  to  be  developed  can  be  tripped  in  the  initial  condition 
with  an  input  to  the  velocity  of  some  kind  instead  of  by  turblence.  The 
Kelvin-Hclmholtz  problem  is  a  prime  example  of  this  where  a  sinusoidal 
velocity  in  the  y  direction  is  used  to  cause  the  instability  to  emerge  in  the 
jet.  The  Euler  equations  are 


dp  i  dpu  |  dpv  |  dpw  n 
dt  ^  dx  '  dy  '  dz 

(50) 

dpu  |  dpu2  +  p  |  dpuv  |  dpuw  n 
dt  dx  '  dy  '  dz 

(51) 

dpv  dpuv  dpv2  +  p  dpwv 
dt  ^  dx  ^  dy  '  dz 

(52) 

dpw  dpuw  dpvw  dpw2  +  p 
dt  '  dx  '  dy  '  dz 

(53) 

dE  du(E  +  p)  dv(E  +  p)  dw(E  +  p) 
dt  +  dx  '  dy  '  dz 

=  0  (54) 

An  assumption  is  made  that  the  air  can  be  represented  as  an  ideal  gas  such 
that  the  specific  heat  leads  to  7  =  Cp/Cv.  This  assumption  also  leads  to  a 
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constant  ideal  gas  constant  R  which  will  be  important  for  normalization.  A 
constant  7  was  chosen  of  1.4.  This  also  causes  a  simple  relationship  between 
the  energy  of  the  problem  and  the  pressure,  density  and  velocities. 


E 


■  P  (  2  1  2  1  2'\ 

+  -  \U  +V  +  W  ) 


(55) 


For  simpler  analysis  the  state  and  fluxes  can  be  stored  into  vectors  so 
that  the  Euler  equation  breaks  down  into  one  vector  equation 


<9u  dF  dG  <9H 

dt  dx  +  dy  dz 


Where  each  vector  represents 


(56) 


p  \ 

(  pll  \ 

(  pV  \ 

(  pw  \ 

pu 

pu 2  +  p 

puv 

puw 

pv 

F  = 

puv 

,G  = 

pv 2  +  p 

■  H  = 

pvw 

pw 

puw 

pvw 

pw2  +  p 

E  ) 

\  u(E  +  p)  ) 

\  v(E  +  p)  ) 

\  w(E  +p)  ) 

(57) 


A  common  choice  for  the  flux  is  the  Lax-Friedrichs  Flux  function  but  this 
is  not  the  best  choice  especially  when  working  in  CFD.  The  AUSM+-up  for 
all  speeds  flux  function  developed  by  Meng-Sing  Lou[14][15]  was  chosen  due 
to  its  capabilities  for  both  high  and  low  speed  flows  and  its  dependence  on 
the  physics  of  the  problem. 


4.1  Normalization 

The  variables  for  this  research  were  all  normalized  by  standard  values.  The 
assumption  of  the  air  being  an  ideal  gas  as  previously  stated  allows  for  the 
constant  specific  heat  and  gas  constants.  A  reference  length  was  chosen 
which  depends  on  the  problem.  For  the  shock  tube  and  Kelvin-Hclmholtz 
instability  the  length  factor  can  be  chosen  indiscriminately  due  to  the  nature 
of  the  problem,  but  for  the  vortex  problem  this  reference  length  plays  great 
weight  in  the  development  of  the  vortex  bursts  initial  condition.  Thus  the 
radius  core  length  rc  was  chosen.  Finally  a  reference  speed  was  chosen  as  the 
speed  of  sound  of  stagnant  air  which  is  a  constant.  This  led  to  the  major 
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variables  being  nondimensionalized  as 


*  x 

X  =  —, 

rc 

*  u 
u  =  — , 

^OO 

_P_ 

Poo 

c*  -  7 

p  7-1’ 


P  = 


y  = 


V  = 


V  = 


e  = 


y 

J 

rc 

v 

®oo 

JL 

Poo' 

e 

a;L 


*  2 

2  =  — , 

rc 

*  w 
w  =  — , 

&00 

T 

rp*  _ 

r  =  — 

To 


(58) 


After  this  point  the  superscript  *  will  be  dropped  and  all  variables  unless 
otherwise  specified  will  be  nondimensional. 


4.2  AUSM+-up  for  all  speeds 

The  AUSM  flux  function  family  was  developed  with  the  goal  of  building  a 
better  flux  function  that  captures  the  accuracy  of  flux  differencing  methods 
with  the  efficiency  of  flux  vector  splitting  methods.  This  method  was  ad¬ 
vanced  substantially  in  [14]  [15]  becoming  first  the  AUSM+  method  and  then 
the  AUSM+-up  method.  The  first  improvement  gave  this  flux  function  the 
ability  to  exactly  capture  a  shock  or  discontinuity  as  well  as  an  improve¬ 
ment  in  accuracy  with  an  easier  adaptation  to  other  conservation  laws.  The 
AUSM+-up  for  all  speeds  method  was  then  created  10  years  later  with  a  goal 
that  as  M  — y  0  the  method  provides  accurate  solutions  therefore  creating  a 
function  that  works  for  all  speeds  improving  the  AUSM+  even  further. 

For  the  development  of  this  method  the  flux  is  broken  up  into  a  convective 
and  pressure  fluxes  such  that 


F  =  mtp  +  P 


(59) 


where 

^  =  (1  ,u,H)T  (60) 


Here  the  left  state  which  corresponds  to  the  state  from  the  left  cell  will  be 
upheld  by  a  subscript  L  and  the  right  state  will  have  a  subscript  R.  The 
mach  number  M  is  found 


Ml/r 


ul/r 

0*1/2 


(61) 
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(62) 


where  ai/2  represents  a  common  speed  of  sound  found  by 

Q-l  +  or 

“1/2  - 

A  common  speed  of  sound  is  very  useful  because  through  the  use  of  a  com¬ 
mon  speed  of  sound  the  exact  capturing  of  both  shocks  and  discontinuities 
becomes  possible. 

The  Mach  number  at  the  interface  is  represented  by  a  function  -M(m\ 
where  m  represents  the  order  of  the  function 


^±(M)  =  t(M±|M|) 

(63) 

yyj  =  i(M  ±  i)2 

(64) 

(4>'  \^)(M)(1T16^||(M)), 

if  |M|  >  1, 
otherwise, 

(65) 

here  f3  represents  a  constant  value.  With  the  Mach  functions  set  the  Mach 
number  at  the  interface  is  then  found  by  applying 


M1/2  —  ^A(Ml)  —  ^#/4x(M  r)  +  Mf 


(66) 


where  Mp  represents  the  pressure  diffusive  terms  effects  on  speed.  The  pres¬ 
sure  diffusive  term  is  represented  by 


K, 


Mp  =  — -max(l  —  crMy  0) 

Ja 


PR~  PL 
Pl/2ay2  ’ 


P 1/2 


PL  +  PR 


(67) 


where  fa  is  a  scaling  function  and  Kp  and  cr  are  constants.  The  mean  local 
rnach  number  (M)  is  found  through 


M2  =  M  +  M«) 

Zal/2 

Which  is  then  used  to  find  a  reference  mach  number 

Mg  =  min  (l,  max  (M2,  M^,))  G  [0, 1] 
which  gives  the  scaling  function  used  above 

/o(M0)  =  M0(2  —  M0)  e  [0, 1] 


(68) 


(69) 
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(70) 


Using  the  new  Mach  number  at  the  interface  and  the  common  speed  of 
sound  the  mass  flow  rate  is  found 


hr  1/2  —  GH/2M1/2 


if  M1/2  >  0, 
otherwise 


(71) 


The  pressure  term  at  the  interface  is  then  found  by  use  of  a  &(n)  function 
which  is  of  varying  degree  n  like  the  Mach  function  previously. 


if  lMl  - 

•/#(^[(±2  —  M)  =1=  IQaM  ^2]  otherwise 


(72) 


where  a  is  a  function  based  on  the  scaling  function  fa.  This  function  is  then 
used  to  develop  the  new  pressure  term  for  the  function  where 


P1/2  =  ^5)(^Il)pl  +  ^J)(M R)pR  -  pu  (73) 

pu  represents  the  velocity  diffusion  term  effecting  low  velocity  flow  and  is 


pu  =  AU^(t)^(5)(PL  +  pR)(faai/2)(uR  -  uL ) 

here  Ku  is  a  user  defined  constant. 

The  previous  variables  a  and  f3  are  then  found  through 


“=n<- 4  +  5^)  €[-14], 


(74) 


(75) 


where  /3  is  left  as  a  constant  and  now  a  has  become  a  function  which  differs 
from  the  AUSM+  method.  Finally  the  flux  is  returned  as 


fl/2 


if /hi/2  >  0 
otherwise  7>i!2 


(76) 


For  supersonic  flows  the  pressure  term  in  the  Mach  number  calculation 
and  the  velocity  term  in  the  pressure  calculation  fall  away  revealing  the 
method  to  be  the  same  as  the  AUSM+.  Therefore  the  method  is  only  being 
adjusted  for  low  speed  flows. 
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4.3  Integration  through  time 

The  DG  method  allows  for  a  superb  method  for  finding  the  derivatives  in 
the  spatial  direction  but  the  problem  remains  with  integrating  these  with 
respect  to  time.  Before  a  scheme  can  be  built  for  this  a  time  step  needs  to  be 
established.  A  reasonable  time  step  is  the  basis  for  any  numerical  calculation. 
Here  the  time  step  can  be  found  by 

(  N2  N4\ 

At  ~  C  (  \max  —  b  IMU00-^"  J  (77) 


where  \max  represents  the  maximum  characteristic  velocity  A max  =  \fypJ~pA 
\/u2  +  v 2  +  w2,  h  is  the  minimum  length,  and  v  is  the  maximum  applicable 
viscosity.  The  maximum  viscosity  coefficient  is  be  found  by  using  a  simple 
relationship  where 


v  = 


A  max 


h 


(78) 


With  the  time  step  established  a  system  for  integrating  through  time 
needs  to  be  found.  The  most  common  and  useful  method  for  doing  this 
is  through  a  Runge-Kutta  integration.  A  total  variation  diminshing  (TVD) 
Runge-Kutta  was  used  in  order  to  remove  any  superfluous  oscillations  caused 
by  the  integration  through  time.  This  is  a  common  method  of  applying  the 
Runge-Kutta  for  a  more  accurate  solution. 

For  this  simulation  a  3rd  order  TVD  Runge-Kutta  integration  was  used 
such  that 


M(d  =  un  +  A tn  ( un ) 
u{2)  =  -un  +  -u(1)  +  -A tn  («(1)) 

4  4  4  v  ’ 

un+1  =  -un  +  -u^  +  -A tU  (u(2)) 

3  3  3  v  ; 


where  1Z  represents  du/dt. 


4.4  Artificial  Viscosity  Sensor 

The  DG  method  has  a  large  advantage  over  traditional  FV  codes  due  to  it 
being  able  to  achieve  a  high  order  of  accuracy.  There  is  a  problem  though. 
Inherent  in  any  high  order  scheme  are  oscillations  created  by  the  Gibb’s 
phenomena.  This  phenomena  can  lead  to  instabilities.  Therefore  a  need  is 
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found  for  a  scheme  to  develop  artificial  viscosity  to  be  used  for  discontinu¬ 
ities.  Traditional  finite  volume  schemes  use  two  different  methods  for  this 
implementation.  The  first  is  to  allow  the  inherent  numerical  dissipation  of 
the  problem  to  remove  these  problems.  This  is  only  applicable  to  linear  first 
order  systems  and  are  not  useful  for  solving  problems.  The  second  is  applied 
to  second  order  and  higher  methods  which  is  accomplished  by  making  the 
scheme  nonlinear.  This  is  done  by  either  applying  a  limiter  to  reduce  the 
order  of  the  simulations  so  that  the  inherent  numerical  dissipation  to  affect 
the  solution  or  by  creating  an  explicit  viscosity  term  based  on  the  state  and 
applying  it.  The  concept  of  limiters  were  dropped  in  order  to  allow  for  a  rep¬ 
resentation  of  the  physics  to  be  the  base  upon  which  the  viscosity  is  applied. 
This  is  artificial  damping. 

A  simpler  way  would  be  use  a  global  viscosity  on  all  cells  much  like  the 
inherent  dissipation  present  in  lower  order  methods  but  this  can  cause  the 
erasing  of  important  small  scale  physics  that  should  not  have  been  affected 
if  this  were  a  real  world  experiment.  If  global  viscosity  is  not  useful,  other 
traditional  codes  apply  a  type  of  limiting  [17]  where  the  DG  method  is  scaled 
back  in  certain  areas  to  lower  the  order  of  accuracy  and  insert  some  of  that 
numerical  dissipation.  This  is  also  not  a  useful  method  because  it  lowers 
the  overall  accuracy  of  the  solution  for  those  cells  even  if  the  cells  are  small 
enough  to  capture  the  ”  sharp”  features.  Therefore  a  sensor  approach  is  much 
more  useful  in  that  it  senses  which  cells  need  to  have  viscosity  applied  for 
them  to  function  without  any  reduction  in  accuracy  for  those  cells.  The 
current  sensor  being  used  was  developed  by  Andreas  Klockner  [12]  as  an 
extension  to  DG  codes  and  can  be  used  with  the  Euler  equations.  Some 
improvements  have  been  made  to  this  code  and  will  be  noted  below. 

The  governing  equations  are  listed  below  where  the  right  hand  side  con¬ 
tains  an  artificial  viscosity  with  u  representing  the  viscosity  coefficient  of 
unknown  value. 

dtp  +  V,  •  (pu)  =  V,  •  (vVxp) ,  (80) 

d t  (pu)  +  V*  •  (u  <8>  (pu))  +  Vxp  =  \7X  ■  (t'V x  (pu)) ,  (81) 

dfE  +  Vo;  •  (u  (A  +  p))  —  Wx  ■  (z'V XE)  (82) 

In  the  development  of  this  sensor  density  was  used  as  the  sensing  variable 
upon  but  it  is  a  choice  of  the  programmer  in  selecting  which  variable  to  sense 
upon.  Here  q  will  represent  the  chosen  state. 

Once  the  choice  is  made  the  Vandermonde  matrix  is  applied  to  transform 
from  a  nodal  basis  into  a  modal  basis.  Here  q  represents  the  given  state  after 
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having  been  transformed  into  a  modal  form.  It  is  assumed  that  the  decay 
can  be  approximated  as 

\qn\  ~  cn~s  (83) 

which  by  taking  the  logarithm 

log|<?n|  ~  log  (c)  -  s  log  (n)  (84) 


gives  an  algebraic  system  of  equations  to  solve.  The  coefficients  are  found  by 
a  least  squares  fit  through  utilization  of  the  normal  equations  where 


AtAx  =  ATb 


(85) 


1  -log(l) 

'loglgil ' 

A  = 

1  — log(2) 

,  b=< 

log  1^1 

1  —  log(n) 

}og\qn\j 

(86) 


s  represents  the  decay  coefficient  and  is  key  in  determining  how  much  vis¬ 
cosity  should  be  applied  to  a  cell.  The  problem  is  this  can  be  misled  very 
easily  by  white  noise  in  the  data  or  for  a  ”kink”  where  certain  modes  fall  to 
zero  leading  to  the  data  being  misinterpreted  and  giving  a  very  high  decay 
coefficient  value. 

To  help  fix  this  two  different  fixes  were  applied  to  the  modes  before  the 
normal  equations  were  solved.  The  first  fix  was  by  applying  skyline  pes- 
simization  which  changes  the  modes  so  that  they  follow  a  monotone  mode 
profile  by  eliminating  spurious  nodes.  Skyline  pessimization  is  implemented 
such  that 


qn  :=  max  |&|  for  n  e  {1,  2, . . . ,  Np  —  1}  (87) 

iG{min(n,ATp— 2),...,ATp— 1} 


Thus  producing  a  new  set  of  modal  coefficients. 

The  small  white  noise  perturbations  are  removed  by  adding  a  represen¬ 
tation  of  scale  to  the  model  in  the  form  of  baseline  decay 


|  K 


r^j 


1  1 


/v-OVp— 1  i  Tl  ^ 

V  2^i=l  ^2 N 


(88) 


Np-1 

£  K\2  =  1 


n= 1 


(89) 
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Which  is  then  added  to  the  original  modal  array  by 

kn|2  :=  \Qn\2  +  IlgjvlllapjSnl2  far  G  {1,  .  .  .  ,  Np  -  1}  (90) 

This  addition  drowns  out  the  floating  point  white  noise  in  the  system  allowing 
for  only  the  data  to  be  examined. 

With  the  value  of  s  determined  it  is  noted  that  an  analogy  can  be  shown 
between  Fourier  decay  and  this  modal  decay  leading  to  the  knowledge  that 
s  ~  1  for  a  discontinuous  solution,  s  =  2  for  a  C°  solution,  s  =  3  for  a  C1 
solution,  etc.  This  leads  to  the  artificial  viscosity  being  determined  by 

(l  se(- 00,1), 

Ks)  =  vq  \  1(1  +  sin(-  (s  -  2)  7t/2))  se[l,3],  (91) 

[o  s  G  (3,  oo). 

where  vq  is  the  maximum  viscosity  coefficient  discussed  previously. 

Their  method  only  included  details  up  to  the  second  dimension  so  in 
order  to  use  Klockner’s  sensor  in  three  dimensional  space  some  new  imple¬ 
mentations  had  to  be  added.  These  changes  caused  both  an  increase  in  the 
speed  of  the  sensor  by  augmenting  the  process  where  the  inverse  Vander¬ 
monde  matrix  was  applied.  A  change  was  implemented  that  moved  the  two 
dimensional  application  from  using  the  Vandermonde  based  off  of  a  two  di¬ 
mensional  system  to  that  based  off  of  a  one  dimensional  system.  The  state 
vector  q  was  split  from  a  vector  of  Np  points  to  a  NfpxNfp  matrix  allowing 
for  application  of  the  Vandermonde  matrix  based  off  of  a  one  dimensional 
system.  This  was  done  through  a  manipulation  of  pointers  in  fortran  90. 
The  norm  was  also  no  longer  used  because  as  the  number  of  dimensions  in¬ 
creased  the  norm  increased  and  with  this  increase  the  norm  could  buffer  out 
some  of  the  important  data  by  making  the  baseline  decay  normalization  too 
big.  Another  problem  exists  because  for  some  data  like,  axial  velocity  in 
the  burst  problem,  the  norm  would  go  to  zero  removing  a  key  piece  of  the 
sensor’s  development.  A  ’’fix”  was  made  to  allow  for  sensing  on  this  data 
and  is  implemented  in  Chapter  4. 

4.5  Improving  the  Efficiency  of  the  Code  to  Increase 
its  Speed 

With  all  these  other  implementations  two  small  but  major  changes  were 
implemented  in  order  to  drastically  increase  the  speed  of  the  code  as  it  was 
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developed  and  moved  from  one  dimension  to  a  three  dimensional  problem. 
These  two  ideas  were  the  creation  of  sparse  matrices  for  stiffness  and  lifting 
matrices  as  well  as  the  implementation  of  parallelization. 

4.5.1  Sparse  Matrix 

The  point  of  sparse  matrix  multiplication  is  to  remove  any  cells  that  are 
zero  or  roughly  equate  to  zero  in  terms  of  double  precision  for  computers. 
These  matrices  had  a  vast  number  of  points  where  the  stiffness  matrix  for  two 
dimensional  problems  had  dimension  [Njp,  Nj]  with  only  Njp  points  having 
values.  Here  Nfp  represents  the  number  of  points  to  a  edge  or  Nfp  =  N  +  1. 
In  terms  of  numbers  this  lays  out  to  a  9th  order  polynomial  having  a  matrix 
of  10,000  points  but  only  1000  points  having  needed  data.  This  data  emerges 
into  patterns  of  numbers  as  seen  in  the  stiffness  and  lifting  matrices  in  the 
figures  below. 


(a)  Sx 


(b)  Sy 


Figure  4:  Location  of  all  non  zero  values  in  stiffness  matrices 
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(cl)  LXa  (b)  Lya 

Figure  5:  Location  of  all  non  zero  values  in  lifting  matrices 


This  is  also  only  for  two  dimenisonal  problems,  the  number  of  wasted 
data  points,  thus  leading  to  pointless  calculations,  only  grows  larger  when 
applied  to  three  dimensional  problems. 

Therefore  a  method  had  to  be  and  was  developed  in  order  to  take  advan¬ 
tage  of  this  and  speed  up  the  code.  A  module  was  created  in  Fortran  that 
searched  the  matrix  for  all  pertinent  data  that  resided  above  a  tolerance  and 
stored  that  into  a  new  data  array.  This  allowed  for  a  speed  up  of  the  code 
due  to  a  smaller  flop  count  in  the  multiplication. 


4.6  Parallel  Processing 

To  obtain  a  faster  code  parallel  processing  was  implemented  such  that  the 
full  power  of  the  processor  could  be  achieved.  Parallel  processing  occurs 
when  a  series  of  loops  is  compressed  and  then  broken  up  and  sent  into  a  user 
defined  number  of  threads  in  a  processor  or  graphics  processing  unit.  For 
these  simulations  a  built  in  method  called  openmp  was  utilized  to  process 
the  parallized  code.  A  simple  experiment  was  done  in  order  to  determine 
the  optimal  number  of  threads  that  should  be  implemented  for  the  xenon 
processor  in  the  workstation  used.  The  data  can  be  seen  in  the  figure  below. 
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Ni=Nj=5 


Ni=Nj=1 5 


(a)  Ni  =  Nj  =  5  cells  (b)  Ni  =  Nj  =  15  cells 

Figure  6:  Wall  time  for  utilizing  multiple  threads  for  the  Isentropic  Vortex 
Test  Case. 

This  data  shows  that  seven  threads  gives  an  optimal  solution  but  this 
could  also  be  weighed  by  the  fact  that  during  these  simulation  the  computer 
had  a  user  logged  into  it.  Nothing  was  running  but  the  eighth  thread  was  oc¬ 
cupied  with  the  operating  system  and  running  the  standard  desktop  functions 
slowing  down  the  system.  Therefore  if  the  user  is  logged  in  seven  threads 
is  optimal  but  for  use  on  something  like  a  cluster  the  maximum  number  of 
threads  or  8  will  be  more  optimal. 


5  Sensor  Improvement 

With  the  substantial  cost  in  terms  of  time  of  the  burst  problem  along  with 
the  different  improvements  made  to  the  artificial  viscosity  sensor  multiple 
test  cases  were  run  for  validation  to  lead  to  the  full  scale  problem.  Previous 
research  accomplished  at  Auburn  University  by  Reitz  [19]  and  Favors  [7]  pro¬ 
vided  insight  into  some  2D  applications.  For  Reitz  it  was  the  implementation 
of  the  9th  order  polynomials  as  well  as  the  importance  of  domain  size  being 
most  favorable  for  speed  and  accuracy  while  Favor’s  thesis  showed  some  of 
the  benefits  that  comes  with  the  Klockner  sensor  used  here  versus  an  entropy 
physics  based  model. 


Auburn  University 
34 


5.1  Test  Cases 


For  the  development  of  the  system  two  separate  two  dimensional  test  cases  as 
well  as  a  three  dimensional  test  case  were  examined  with  the  sensor.  These 
cases  are  a  shock  tube  adjusted  to  two  dimensional  space  and  the  Kelvin- 
Helmholtz  phenomena.  Then  before  the  implementation  of  the  full  three 
dimensional  problem  the  shock  tube  was  extrapolated  to  three  dimensions  to 
validate  the  changes  that  were  made  to  the  sensor  as  well  as  work  with  some 
new  changes  specifically  in  how  the  baseline  decay  is  applied. 

5.1.1  Kelvin-Helmholtz 

The  Kelvin  Helmholtz  phenomena  was  used  as  a  simple  2D  model  in  order 
to  validate  the  changes  as  well  as  show  an  increase  in  the  speed  of  the  new 
implementation  of  the  sensor.  Another  advantage  to  using  this  as  a  check,  is 
that  the  Kelvin-Helmholtz  phenomena  is  a  naturally  occurring  phenomena 
that  can  occur  in  wakes  giving  one  more  example  of  the  DG  method  modeling 
wake  flow  components.  This  test  case  also  gives  a  special  situation  where 
the  sensor  can  be  checked  against  numerical  instabilities  that  are  developed 
through  time  and  are  not  abrupt  like  a  shock  is. 

In  creating  the  initial  state  U  constant  values  are  assumed  for  p  =  1  and 
p  =  l/y.  Constants  e  =  vr/15,  5  =  0. 1/4.9,  and  M  =  0.25  are  used  in  the 
development  of  the  velocity  profile  for  the  system  which  is  found  through 

tanh  /37rA~yj  ,  y  >  tt 
tanh  j  ;  else 

u  =  0.5  (s  +  1)  M 

v  =  M<5  (1.  sin  (1.  x/3.)  —  1.  sin  (2.  x/3.)  +  1.  sin  (3.  x/3.)) 

The  vertical  velocity  here  is  developed  with  the  purpose  of  perturbing  the 
first  three  modes  to  allow  for  consistently  reproducible  simulations  in  order 
to  compare  results  back  and  forth.  A  fifth  state  is  added  to  the  governing 
equation  concerning  the  concentration  factor  s.  This  state  is  used  both  in 
the  sensor  and  in  the  plotting  of  the  initial  conditions  as  seen  in  Figure  7. 
This  allows  for  an  easy  comparison  to  other  codes  and  experiments.  The 
bounds  for  this  problem  are  x  G  [0, 47r]  and  y  G  [ — 7T,  37t]  . 
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Figure  7:  Initial  condition  showing  the  concentration  s 

5.1.2  Sod  Shock  Tube 

The  burst  problem  also  deals  with  a  direct  shock  like  discontinuity  when  the 
vortices  core  radius  changes  due  to  the  burst.  Therefore  an  examination  of  a 
shock  tube  provides  ample  benefits  in  how  the  sensor  handles  such  an  abrupt 
change  in  the  state.  A  shock  tube  is  established  in  experiments  by  creating 
a  cylinder  with  two  different  states  for  pressure  and  density  on  either  side 
of  a  burst  disk.  At  time  t=0  the  burst  disk  is  then  popped  and  the  gases 
mix  sending  a  shock  wave  and  other  cascading  effects  through  the  system. 
This  problem  from  a  CFD  sense  is  a  simple  one  dimensional  problem  that 
can  easily  be  extrapolated  to  three  dimensional  space. 

The  problem  here  examines  the  common  test  case  of  the  Sod  shock  tube 
and  is  set  up  for  multiple  dimensions  by  choosing  a  single  direction  and  then 
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declaring  the  location  of  the  burst  disk  between  the  two  ends  of  the  domain 
to  exist  in  this  direction.  For  this  simulation  a  bounds  of  [—1, 1]  for  both  x 
and  y.  The  computational  burst  disk  was  placed  at  0  which  is  the  center. 
The  values  on  either  side  are 


(  1.  x  <  0. 
p  ~  \  0.125  x  >  0. 

f  1.  x  <  0. 

(  0.1  x  >  0. 

and  an  image  of  the  initial  condition  can  be  found  in  Figure  8. 
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Figure  8:  Initial  Condition  of  density  (p) 

With  the  sensor  fully  validated  the  shock  tube  was  then  expanded  into 
three  dimensions  in  order  to  troubleshoot  as  well  as  develop  the  extrapolation 
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of  the  sensor  into  three  dimensional  space.  The  same  format  for  the  develop¬ 
ment  of  the  initial  condition  was  followed  for  converting  the  two  dimensional 
shock  tribe  to  a  three  dimensional  shock  tube. 

5.2  Results  for  Sensor  Validation 

5.2.1  Sensor  Validation  in  Two  Dimensions 

The  modal  sensor  is  a  very  useful  tool,  but  for  an  optimal  use  in  three  di¬ 
mensions  it  needs  to  be  retooled  to  increase  it’s  speed.  This  is  an  important 
concept  when  running  three  dimensional  simulations  which  are  already  very 
costly  concerning  time.  Before  moving  into  the  three  dimensional  problems, 
the  sensor  was  first  altered  in  the  two  dimensional  space  in  order  to  trou¬ 
bleshoot  as  well  as  see  if  it  altered  the  sensor’s  abilities.  These  changes  can 
be  found  in  algorithms  1  and  2.  Here  the  original  algorithm,  algorithm  1, 
has  a  matrix  multiplication  that  costs  Njp  flops.  This  is  due  to  it  involving 
a  matrix  multiplication  of  a  NjpxNjp  matrix  to  a  Nfpxl.  The  method  is 
altered  by  using  the  inverse  of  the  Vandermonde  matrix  based  off  of  a  one 
dimensional  system  instead  of  the  current  two  dimension  system.  This  Van¬ 
dermonde  matrix  was  then  applied  by  matrix  multiplication  in  the  x  or  y 
direction  while  looping  through  all  the  nodes.  This  resulted  in  2  NfpxNfp  by 
Nfpx  1  calculations  inside  of  a  loop  over  Nfp  creating  a  flop  count  of  2  *  Njp 
for  the  two  dimensional  case.  For  the  3D  case  the  savings  are  bolstered  even 
higher.  This  case  used  to  cost  flops  on  the  order  of  N®p  but  with  the  use  of 
the  smaller  Vandermonde  matrix  it  only  costs  3  *  Nj  flops.  This  translates 
to  a  savings  in  two  dimensional  of  Nfp/2  and  savings  in  three  dimensional 
space  of  Nj/ 3  flops  which  shows  exponential  savings  between  the  two  cases. 

With  proof  that  the  sensor  provides  a  substantial  increase  in  the  speed 
of  the  simulation  the  validity  of  the  new  application  of  the  modal  sensor 
was  checked  using  the  Kelvin-Helmholtz  problem.  This  was  done  to  prove 
the  modes  were  still  the  same.  The  kinetic  energy  in  the  Kelvin-Helmholtz 
instability  is  a  constant.  Therefore  the  kinetic  energy  was  integrated  using 
the  gauss-lobatto  quadrature  and  then  was  normalized  by  the  initial  kinetic 
energy.  This  was  plotted  in  Figure  9.  Global  viscosity,  or  the  application  of 
viscosity  to  all  locations,  was  also  plotted  to  show  how  much  more  effective 
using  a  sensor  is.  The  figure  shows  both  the  new  approach  and  old  approach 
giving  an  identical  integrated  kinetic  energy  over  time  while  providing  leaps 
and  bounds  better  results  than  a  global  viscosity  method  proving  the  validity 
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of  the  sensor  while  also  showing  how  it  obtains  similar  results. 

The  sensor  was  then  compared  as  it  activated  in  both  x  and  y  direction 
as  seen  in  Figures  10  through  13.  As  can  be  seen  in  the  figures  the  sensor 
turned  on  in  the  same  spots  with  the  exact  same  weight  for  both  the  new  and 
old  implementations  which  aligns  with  information  found  in  Figure  9.  This 
all  logically  follows  since  the  modes  are  the  same  regardless  of  the  method 
for  finding  them. 


Figure  9:  Total  kinetic  energy  over  time 
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(a)  ex  New  (b)  ex  Old 

Figure  10:  Sensor  in  X  Direction  at  time  49.49  s 


5  10  15  20  25  5  10  15  20  25 


(a)  ey  New  (b)  ey  Old 

Figure  11:  Sensor  in  Y  Direction  at  time  49.49  s 
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(a)  ex  New  (b)  ex  Old 

Figure  12:  Sensor  in  X  Direction  at  time  76.99  s 


(a)  ev  New  (b)  ey  Old 

Figure  13:  Sensor  in  Y  Direction  at  time  76.99  s 
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(a)  Kelvin-Helmholtz  New  (b)  Kelvin  Helmholtz  Old 

Figure  14:  Comparison  of  the  two  states  at  time  76.99  s 

With  validation  that  the  new  implementation  returns  the  same  amount 
of  artificial  viscosity  for  both  implementations  a  cost  comparison  was  done 
by  utilizing  both  the  Sod  shock  tube  and  the  Kelvin-Helmholtz  phenomena. 
Tables  1-3  show  the  new  implementation  of  the  sensor  saves  time  for  each 
case.  In  each  of  these  tables  the  total  time  for  running  each  simulation  is 
what  is  recorded.  Tables  2  and  3  also  compare  an  increase  in  the  order  of 
the  polynomial  and  how  that  affects  the  new  method.  For  each  fo  the  sod 
shock  tube  test  cases  a  constant  At  was  used  in  order  to  compare  the  results 
for  the  varying  polynomial  order.  For  the  data  stored  in  Table  2  the  At  was 
0.00005  and  for  Table  3  it  was  0.0002.  For  the  Kelvin-Helmholtz  there  is 
a  4.0858%  savings  in  time  while  for  the  Sod  shock  tube  there  is  a  3.9611% 
savings  in  wall  clock  time  for  a  ninth  order  polynomial.  This  shows  a  decrease 
in  time  but  questions  remain  due  to  the  small  difference  when  comparing  the 
sensor  for  the  test  cases  of  table  2.  Based  on  polynomial  order  the  decrease 
in  time  varies  from  5.318  s.  Therefore  with  the  domain  specifications  of 
Ni  =  99  Nj  =  1  the  order  of  the  polynomial  was  increased  in  order  to  judge 
whether  or  not  the  method  did  allow  for  more  savings  as  polynomial  order 
increased.  For  the  tenth  order  polynomial  the  savings  amounted  to  7.278  s  or 
5.10  %,  the  eleventh  order  polynomial  is  11.015  s  or  6.24  %,  the  twelfth  order 
polynomial  is  15.873  s  or  7.83  %,  and  the  thirteenth  order  polynomial  gives 
22.17  s  or  9.06  %  in  savings.  This  proves  that  as  the  order  of  the  polynomial 
is  increased  the  savings  also  increase.  Since  the  sensor  is  implemented  for 
each  cell  and  the  new  implementation  saves  on  the  order  of  Nfp/ 2  flops  per 
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iteration  the  number  of  cells  was  increased  to  Nt  =  50  Nj  =  50  in  order  to 
see  how  substantial  the  savings  can  be  for  more  realistic  two  dimensional 
meshes. 

The  results  for  these  new  runs  can  be  seen  in  Table  3.  Here  there  is  a 
savings  in  time  of  14.256  s  or  2.20%  for  the  ninth  order  polynomial,  30.507  s 
or  3.79  %  for  the  tenth  order  polynomial,  41.649  s  or  4.09  %  for  the  eleventh 
order  polynomial,  68.24  s  or  5.58  %for  the  twelfth  order  polynomial,  and 
96.433  s  or  6.31  %  for  the  thirteenth  order  polynomial.  This  proves  the 
trend  found  in  the  previous  shock  tube  experiment  where  as  the  order  of 
the  polynomial  is  increased  there  are  more  savings  from  the  augmentation  to 
the  sensor.  The  decrease  in  the  percentage  of  time  can  be  attributed  to  the 
increase  in  mesh  size  and  the  sensor  not  being  the  most  costly  part  of  the 
simulations. 


Algorithm  1  Old  Sensor  Implementation 

1:  for  i  —  1,  Ni  do 

2:  for  j  —  1,  Nj  do 

3:  Q  =  U(:,5,i,j)/U(:,l,i,j) 

4:  Q  —  abs ( rnatmul (inn V,  Q(:,  1))  +  l.Oe  —  17)  >  Njp  flops 

5:  for  in  =  1,  Nfp  do 

6:  qx(iii)  =  sqrt(sum {Q(KceU(\,  in),  1)  *  *2.)) 

7:  qy(iii)  =  sqrt(sum (Q(Kceu(iii, :),  1)  *  *2.)) 

8:  end  for 

9:  : 

10:  end  for 

11:  end  for 


Auburn  University 
43 


Algorithm  2  New  Sensor  Implementation 
1:  for  i  —  1,  iVj  do 
2:  for  j  =  1,  Nj  do 

3:  call  wrapper  (din,  U,  uhat ) 

4:  =  apply-along-all(<im,  invVID,  uhat)  >  2  *  lV|p  Flops 

5:  um  =  nm  *  "Um 

6:  gxl  =  sqrt(sum(wm,  2)) 

7:  gyl  =  sqrt(surn(iim,  1)) 

8:  : 

9:  end  for 

10:  end  for 


Sensor  Implementation 

Wall  Clock  Time 

New  Method 

545.3109  s 

Old  Method 

568.5409  s 

Table  1:  Kelvin  Helmholtz,  N=9,  Ni=27,  Nj=18 


Wall  Clock  Time 

Sensor 

N=9 

N=10 

N=ll 

N=12 

N=13 

New  Method 

107.232  s 

135.144  s 

165.39  s 

204.752  s 

222.661  s 

Old  Method 

112.55  s 

142.422  s 

176.405  s 

220.393  s 

243.831  s 

Table  2:  Sod  Shock  Tube,  Ni=99,  Nj=l 


Wall  Clock  Time 

Sensor 

N=9 

N=10 

N=ll 

N=12 

N=13 

New  Method 

634.739  s 

775.489  s 

977.014  s 

1154.97  s 

1431.889  s 

Old  Method 

648.995  s 

805.996  s 

1018.663  s 

1223.21  s 

1528.322  s 

Table  3:  Sod  Shock  Tube,  Ni=50,  Nj=50 
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5.2.2  Implementation  into  Three  Dimensions 

With  both  the  validity  of  the  new  approach  to  the  modal  sensor  proven  as 
well  as  the  savings  shown  to  be  substantial,  the  next  step  was  to  extrapolate 
the  sensor  into  three  dimensional  space.  When  applying  the  sensor  to  the 
vortex  bursting  problem  the  sensor  repeatedly  failed  due  to  it  not  activating 
fast  enough  on  the  data.  Both  density  and  energy  were  sensed  to  no  avail. 
Therefore  the  sensor  was  further  augmented  to  look  for  the  highest  value  over 
all  the  states  and  apply  that  to  the  relative  cell  in  all  3  directions.  A  problem 
emerges  when  applying  this  method  though  due  to  the  norm  for  the  state  pw 
will  be  zero,  therefore  the  sensor  may  turn  on  due  to  the  inherent  floating 
point  noise.  This  can  ”  wash  out”  important  aspects  fo  the  physics  of  the 
problem.  This  meant  that  a  new  method  for  finding  the  stabilizing  factor 
was  needed  for  the  application  of  the  baseline  decay  to  make  the  sensor  more 
robust  on  all  states  even  if  it  removed  some  of  the  simplicity  of  the  sensor. 
To  validate  this  new  approach  the  Sod  shock  tube  was  returned  to  due  to  its 
simplicity  as  well  as  the  speed  with  which  the  problem  can  be  run. 

The  exact,  analytical  answers  for  the  Sod  shock  tube  problem  were  found 
from  a  NASA  website[20]  and  used  to  tabulate  RMS  errors  for  each  simu¬ 
lation.  The  new  scale  for  the  norm  was  originally  developed  as  something 
manually  entered  by  the  user  for  the  simulations.  For  these  simulations 
represented  in  Figures  15-17  the  norm  for  each  state  can  be  found  in  the 
following  vector  ynorm  =  [3.0,0.005,0.005,0.005,5.0].  The  values  for  states 
two  through  four  were  chosen  in  order  to  see  the  level  of  impact  the  velocity 
would  have.  The  values  for  the  first  and  fifth  state  were  chosen  to  be  close 
to  their  norm  values  for  a  one  dimensional  case. 

Returning  to  the  problem  the  test  cases  were  run  with  an  activation  just 
on  density  and  then  on  all  states  in  order  to  provide  a  good  comparison  of 
the  data.  Table  4  shows  for  both  cases  an  RMS  error  of  less  than  0.02  was 
found.  Figure  18  shows  that  by  sensing  on  the  whole  state  the  results  are 
increasingly  smoothed  over  just  sensing  on  one  state  due  to  the  increase  in 
the  activation  of  the  sensor.  This  can  be  further  seen  from  Figures  15-17 
which  depicts  a  comparison  of  the  sensor  being  turned  on  through  time  as 
the  shock  is  propagated  in  any  of  the  three  directions.  Image  (b)  in  each  of 
the  Figures  shows  that  the  sensor  actively  tracks  just  the  shock  and  keeps 
a  good  marker  on  it  when  just  sensing  on  density  while  Image  (a)  in  the 
Figures  shows  a  more  liberal  application  of  the  viscosity  over  a  wider  range 
of  data. 
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Values  Sensed  On 

RMS  Error 

Sense  on  p 

0.0078 

Sense  on  all  states 

0.0142 

Table  4:  RMS  Error  for  Sod  Shock  Tube  3D 
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Cells  in  x  Direction 

(a)  X  Direction  All 
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Cells  in  x  Direction 


(b)  X  Direction  Rho 


Figure  15:  Sensor  through  time  X 
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Cells  in  y  Direction 


(a)  Y  Direction  All 


(b)  Y  Direction  Rho 


Figure  16:  Sensor  through  time  Y 
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(a)  Z  Direction  All  (b)  Z  Direction  Rho 

Figure  17:  Sensor  through  time  Z 
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Figure  18:  Comparison  of  Sensor  using  all  states  or  one 
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Figure  19:  Zoomed  in  comparison  of  the  baseline  decay  variation 

As  previously  stated  these  simulations  were  all  done  with  a  predetermined 
scale  of  ynorm  =  [3.0,0.005,0.005,0.005,5.0].  This  raised  a  question  though 
of  whether  or  not  this  was  the  best  choice.  The  simulation  was  then  rerun 
with  a  new  scale  of  ynorm  =  [3.0,  0.5,  0.5,  0.5,  5.0].  A  comparison  of  the  shock 
tube  sensing  on  all  five  state  can  be  seen  in  Figure  20.  These  images  show 
the  sensor  being  activated  much  more  sparsely  due  to  the  higher  value  on 
the  scale  for  the  velocity  states.  This  follows  the  logic  of  the  scale  being  used 
to  balance  the  baseline  decay  to  remove  floating  point  noise.  The  lower  the 
value  for  the  norm  the  higher  chance  there  is  for  noise  to  affect  the  sensor. 
Table  5  validates  the  previous  reasoning  by  showing  a  much  smaller  RMS 
value  for  the  new  implementation  of  the  scale  versus  the  old  one  and  being 
much  more  in  line  with  RMS  value  for  the  sensor  only  utilizing  density. 
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Figure  20:  Comparison  of  different  preselected  weighting  vectors 


Values  Sensed  On 

RMS  Error 

Velocity  norm  0.005 

0.0142 

Velocity  norm  0.5 

0.0088 

Table  5:  RMS  Error  for  Sod  Shock  Tube  3D 

All  this  bears  the  question  of  how  to  derive  a  system  to  accurately  pre¬ 
determine  values  for  this  new  scale.  As  shown  when  discussing  the  sensor  a 
value  for  s  of  1  or  lower  leads  to  the  maximum  application  of  artificial  viscos¬ 
ity.  Whereas  a  value  for  s  of  3  lead  to  the  least  amount  of  artificial  viscosity 
applied.  Any  values  of  s  over  3  causes  no  artificial  viscosity  to  be  applied. 
The  question  is  then  posed  that  for  basic  phenomena  such  as  floating  point 
noise  or  a  step  function  what  is  the  maximum  value  that  can  be  used  to  trip 
the  sensor  for  either  of  these  situations. 

The  previous  method  for  deriving  the  scale  was  also  found  to  cause  dif¬ 
ferent  values  for  a  step  function  based  at  0  and  a  step  function  based  at  any 
other  height.  In  order  to  incorporate  this  the  step  function  to  be  used  will 
be  based  at  one.  In  order  to  determine  this  the  Erst  thing  that  was  done 
was  to  choose  a  set  of  s  values  to  examine.  For  this  research  the  values  of 
s  =  [1,  2,  3]  were  chosen.  Then  three  separate  polynomial  orders  were  chosen 
to  be  used  as  a  framework  for  the  solution.  The  polynomial  orders  chosen 
were  x  =  [5,  9, 15].  Initial  guesses  were  found  for  these  polynomial  orders  in 
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order  to  obtain  the  desired  s  value.  These  values  can  be  seen  in  Equation  96 
where  the  rows  correspond  to  the  distinct  s  values  and  the  columns  to  the 
polynomial  orders. 


0.0034  0.012  0.066 
0.019  0.22  6.6 

0.071  1.72  180 


(96) 


With  the  initial  conditions  built  a  MATLAB  script  was  run  that  created 
an  initial  guess  by  utilizing  this  data  and  an  intrinsic  interpolation  function  in 
MATLAB.  Then  a  minimizing  function  was  used  called  FZERO.  The  function 
being  fed  to  FZERO  would  take  the  polynomial  order  and  determine  an  initial 
condition  based  off  of  the  desired  test  case  and  then  run  the  sensor  on  it. 
With  the  value  of  s  obtained  the  function  would  then  determine  the  cost 
function  with  s desired,  ~  s sensor-  The  system  would  continually  run  until  it 
achieved  an  s  value  within  a  tolerance  to  the  desired  value.  These  values  for 
both  the  step  function  as  well  as  noise  can  be  seen  below  in  Figures  21  and  22. 
In  these  figure’s  the  average  value  needed  to  activate  the  sensor  was  divided 
by  the  amplitude  to  allow  for  an  easy  application  to  different  problems. 
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Figure  21:  The  minimum  necessary  normalized  average  to  trip  a  specific  s 
value  for  noise 
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Figure  22:  The  minimum  necessary  normalized  average  to  trip  a  specific  s 
value  for  a  step  function 

As  previously  stated  a  standard  polynomial  order  of  9  was  used.  The 
values  for  a  9th  order  polynomial  fall  in  Tables  6.  In  order  to  prevent  the 
sensor  from  turning  on  too  much  an  average  of  the  values  for  s  =  2  for  Noise 
and  s  =  3  for  the  step  were  used.  This  average  could  then  by  multiplied  by 
the  amplitude  of  the  desired  state  in  order  . 
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s 

Noise 

Step 

1 

11.844 

0.83421 

2 

217.75 

15.425 

3 

1717.9 

121.47 

Table  6:  s  values  for  9th  order  polynomials 

The  shock  tube  was  then  rerun  with  these  new  scales  to  determine  their 
accuracy.  Figure  23  and  Table  7  show  that  the  new  RMS  error  as  well  as  a 
plot  of  the  new  density  are  nearly  identical  to  the  values  previously  decided 
upon.  These  show  that  the  new  scale  function  will  provide  a  good  result 
while  removing  any  guess  work  previously  being  utilized.  This  is  a  huge 
improvement  in  that  the  sensor  is  more  applicable  to  any  type  of  state  at  a 
relatively  small  cost  of  being  now  based  on  the  order  of  the  polynomial  while 
still  retaining  it’s  independence  of  the  problem  which  many  other  sensors 
struggle  with. 
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Figure  23:  Using  the  new  scale  derivation  to  replace  yn0rm 


Values  Sensed  On 

RMS  Error 

Utilizing  built  in  scale 

0.0093 

Velocity  norm  0.5 

0.0088 

Table  7:  RMS  Error  for  Sod  Shock  Tube  3D 
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6  Three  Dimensional  Vortex  Problems 


6.1  Generation  of  the  Initial  Condition 

The  burst  vortex  is  a  complicated  problem  when  generating  the  initial  con¬ 
dition.  It  requires  pseudo  steady  state  simulations  as  well  as  extrusions  for 
a  two  dimensional  vortex  into  a  three  dimensional  field. 

6.2  Velocity  Profile 

In  order  to  establish  the  real  world  application  of  the  bursting  an  accurate 
representation  of  the  velocity  profile  needs  to  be  made.  For  this  problem  a 
Gaussian  vortex  velocity  profile  was  chosen  that  is  developed  as 

^  )  (97) 

where  Figure  24  shows  the  non-dimensional  velocity  profiles  for  the  two  vor¬ 
tices  with  cores  sizes  of  rCl  =  2.6  and  rC2  =  5.2. 
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Figure  24:  vq  for  the  vortices  of  two  different  core  sizes 


With  the  velocity  profile  built  in  radial  coordinates  it  then  had  to  be 
shifted  to  Cartesian  coordinates.  This  was  done  by  first  transforming  the 
velocity  equation  into  a  velocity  potential  equation  by 


/  = 


ve 

r 


(98) 


There  is  a  problem  with  the  profile  though.  At  the  center  r  =  0  the 
current  velocity  profile  would  reach  a  singularity,  therefore  a  high  order  alge¬ 
braic  model  was  used  only  at  r  —  0.  The  velocity  potential  for  all  radii  can 
then  be  found  from 


Ia 

2  7ra 


(99) 
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(r/re)2+ 27 


f  —  j  rci  V(^/rc^2+7) 

,c  ,2  fl.O  -e~^r^ 

rCi (r/rci )  V 


(100) 


where  T0  represents  the  circulation. 


6.3  Boundary  Condition 

With  a  velocity  profile  the  boundary  conditions  then  had  to  be  developed. 
This  way  the  pseduo-steady  state  simulations  could  be  run  for  the  two  di¬ 
mensional  slices.  This  was  done  with  symmetric  boundary  conditions  in  the 
transverse  directions  [22],  These  are  applied  by  imagining  a  ghost  set  of  do¬ 
mains  surrounding  the  vortex  and  each  ’’ghost”  domain  bears  a  reflection  of 
the  vortex.  This  reflection  is  captured  in  rings  around  the  current  domain. 
Each  ring  contains  more  reflections  of  the  domain  due  to  a  larger  radius  away 
from  the  original  domain  but,  the  further  out  the  ring  is  the  less  the  velocity 
of  this  ’’ghost”  domain  influences  the  development  of  the  velocity  profile. 
The  symmetric  boundary  condition  allows  no  velocity  to  cross  its  boundaries 
therefore  the  domain  of  the  problem  is  of  the  utmost  importance. 
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Figure  25:  Example  of  Symmetric  Boundary  Condition  on  a  Vortex 

For  this  problem  25  rings  were  generated  and  then  a  Van  Wijngaarden 
transformation  was  used  in  order  to  accelerate  the  series  and  find  a  more  exact 
solution  for  the  velocity  than  simply  doing  a  summation.  This  is  performed 
by  first  performing  an  Euler’s  transform  which  is  accomplished  by  doing  a 
summation  of  the  array  in  question 


k 

=  ^(-lfan  (101) 

n= 0 

here  s  is  a  matrix  of  data  for  computing  the  Euler’s  transformation.  The 
velocity  summations  are  then  stored  in  the  top  row  of  data  of  s  at  which 
point  the  Euler’s  transformation  is  performed  where 

_  sj,k  +  sj,k+ 1 
sj+l,k  ~  2 

Van  Wijngaarden’s  contribution  lay  in  noticing  that  it  was  pointless  to  takes 
this  all  the  way  to  the  end  of  the  data.  Rather  by  stopping  two-thirds 
of  the  way  through  the  transformation  a  more  accurate  value  for  the  true 
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summation  was  found.  This  contribution  allowed  for  less  rings  to  be  used 
while  obtaining  much  more  accurate  results  for  the  velocity. 


6.4  Two  Dimensional  Initial  Condition 

While  the  velocity  field  is  well  represented  and  the  boundary  conditions  have 
been  derived,  there  is  no  simple  way  to  derive  the  pressure  and  density  from 
the  velocity.  These  values  are  intrinsically  linked  in  the  Euler  equations.  This 
led  to  a  need  to  develop  a  pseudo-steady  two  dimensional  initialization  algo¬ 
rithm.  This  is  accomplished  by  creating  a  standard  flow  field  with  uniform 
density  and  pressure  with  no  velocity  and  then  slowly  introducing  the  true 
velocity  field  over  a  period  of  time  and  then  running  the  simulation  till  the 
residual  reaches  three  orders  of  magnitude  below  it’s  peak  to  show  a  steady 
state  solution. 

The  introduction  of  the  velocity  is  done  through  a  factor  based  on  a 
cosine  function  going  from  zero  to  one.  After  a  certain  set  time  that  factor 
stays  at  one  so  that  the  initial  condition  can  come  to  a  constant  value.  This 
introduction  equation  is 

ft  7 T  7T\ 

s  =  cos  (  - * -  (103) 

v  100  2  2j  v  ’ 


so  at  a  time  of  one  hundred  seconds  the  value  for  s  remains  at  one. 

To  achieve  the  introduction  of  the  velocity  slowly  into  the  right  hand  side 
of  the  solution  is  as  follows. 


rhs  (:,  2)  =  rhs  2)  +  ^  *  yU  (:,  1  ,i,j)  *  (^ut  (:,M,j) 

rhs  (:,  3)  =  rhs  3)  +  ^  *  (u  (:,1  ,i,j)  *  (vt(:,l ,i,j) 

rhs  (:,  4)  =  rhs  4)  +  ^  *  (u  (:,2 ,i,j)  *  (ut(:,l,i,j) 

+  ^  *  (u  (:,  3,  i,j)  *  (vt(:7l,i,j) 


U( 

A  2,i,  j) 

U( 

),l,bj) 
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3,2  ,i,j) 

U( 

3,1  ,i,j) 

U( 

3,3  ,i,j) 

U( 

3,1  u,i) 

(106) 

Here  i  and  j  corresponds  to  the  cells  in  the  x  and  y  direction,  r  is  a  constant 
held  to  one  half  here  so  that  there  is  more  weight  on  the  velocity  inclusion, 
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and  ut  and  vt  represent  the  true  velocities  in  the  x  and  y  direction  respectively 

This 


6.5  Three  Dimensional  Extrusion 

The  two  two  dimensional  vortices  were  then  generated  with  the  differing  core 
radii  previously  mentioned.  The  ratio  of  the  two  core  radii  was  altered  in 
the  exponential  term  of  Equation  100.  This  way  the  strength  of  the  vortex 
was  unaffected  by  the  changing  core  radius  allowing  for  the  circulation  to 
be  accurately  maintained  and  the  vortices  to  be  merged.  With  the  two 
dimensional  initial  conditions  made  all  that  was  left  was  the  generation  of 
the  three  dimensional  model. 

The  two  dimensional  slices  were  merged  by  using  a  sinusoid  function 
going  from  [1,-1]  for  5.  <  z.  <  10.  and  [—1,1]  for  45.  <  z  <  50..  This  is 
represented  by  the  function  theta  and  can  be  seen  below. 
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(107) 

(108) 


This  generates  the  initial  condition.  The  velocity  in  the  transverse  di¬ 
rections  was  then  perturbed  using  a  random  perturbation  applied  as  u  — 
Wbase  (1  +  eh)  where  e  =  lxlO-2  and  ||u||  <  1.  This  trips  the  helical  insta¬ 
bility  in  the  problem.  The  boundary  condition  for  the  axial  direction  of  the 
vortex  bursting  problem  was  also  a  symmetric  boundary  condition  which  due 
to  the  velocity  being  independent  of  z  is  more  like  a  reflecting  wall  boundary 
condition. 


6.6  Results 

The  CFL  number  plays  a  huge  role  in  the  development  of  the  time  step  for 
all  CFD  applications.  For  the  two  dimensional  cases  this  CFL  number  can 
approach  one  for  most  problems  but  this  is  not  the  case  for  three  dimensional 
problems.  In  order  to  determine  the  appropriate  CFL  number,  the  number 
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was  varied  and  the  three  separate  test  cases  were  run.  These  test  cases  had 
varying  CFL  numbers  of  CFL  =  [0.25,  0.3,  0.35]  with  a  grid  size  of  20x20x20. 
The  simulation  corresponding  to  the  CFL  number  of  0.35  failed,  therefore  to 
keep  the  code  running  as  fast  as  possible  the  CFL  number  was  established 
at  0.3. 

The  minimum  pressure  is  plotted  for  varying  grid  sizes  and  polynomial 
order  in  Figures  26-31.  The  plots  show  that  the  burst  follows  the  same 
path  as  the  vortex  bursting  pressure  held  found  in  Figure  32  [16].  These 
figures  all  show  a  constant  velocity  in  the  pressure  wave  collapsing  into  each 
other.  In  order  to  compare  the  time  the  normalization  was  shifted  in  post 
processing  by  a  conversion  factor  of  vgmax/(2an).  The  burst  is  identified  by 
the  sharp  increase  in  pressure  as  seen  in  each  of  the  Figures  and  is  found 
at  a  constant  location.  This  lines  up  with  the  results  in  the  Moet  et.  al. 
paper  where  they  found  that  a  constant  propagation  speed  existed  for  all 
the  test  cases  and  the  burst  occurs  at  an  exact  time  based  off  the  ratio 
of  the  core  radii.  There  is  a  difference  in  the  amplitude  of  the  pressure 
waves  but  this  is  understandable  due  to  the  difference  in  the  normalization 
techniques.  This  project  was  normalized  by  utilizing  the  dynamic  pressure 
as  the  normalization  factor  where  the  Moet  paper  uses. 

P  Pmin 

p  =  - 

Poo  Pmin 

The  grid  size  and  polynomial  order  where  varied  in  order  to  find  the 
most  accurate  solution  and  see  how  each  change  impacts  the  solution.  The 
Figures  show  that  the  amount  of  cells  in  the  axial  direction  do  not  have  a  large 
influence  on  the  pressure  wave.  Instead  the  transverse  cells  hold  the  most 
influence.  This  is  logical  because  since  this  plot  is  based  on  the  minimum 
pressure  and  plotted  along  the  axial  direction.  A  very  interesting  note  is  that 
even  though  the  13th  order  polynomial  used  less  cells  and  degrees  of  freedom 
it  obtained  an  equally  accurate  version  of  the  solution  when  compared  to  the 
paper  by  Moet  et.  al.  The  downside  to  the  this  implementation  was  in  the 
time  needed  to  run  the  solution.  Even  though  it  had  the  least  degree’s  of 
freedom  the  solution  took  the  longest  to  run. 


(109) 
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Figure  26:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  9th  order  with  a  grid  Ni=Nj= 20  and  N^=2h 
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Figure  27:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  9th  order  with  a  grid  Ni—Nj— 25  and  Nk= 30 


Figure  28:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  9th  order  with  a  grid  Ni—Nj— 25  and  Nk= 35 
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Figure  29:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  9th  order  with  a  grid  Nt—Nj— 25  and  Nk= 40 


Figure  30:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  9th  order  with  a  grid  Nt—Nj— 25  and  Nk= 45 
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Figure  31:  Plot  of  the  minimum  pressure  evolving  through  time  for  a  poly¬ 
nomial  of  the  13th  order  with  a  grid  size  of  Ni—Nj— 17  and  Nk= 25 


Figure  32:  Temporal  evolution  of  the  profile  of  minimum  pressure  in  the 
vortex  core  for  simulation  DNS2[16] 

Since  the  pressure  plot  had  poor  quality  when  the  transverse  number  of 
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cells  was  below  25  in  each  direction  the  vorticity  for  those  solutions  was  not 
plotted.  The  vorticity  magnitude  is  then  plotted  for  each  of  the  remaining 
cases  in  Figures  33-37.  They  are  listed  with  their  time  being  changed  to  the 
scale  of  the  paper  to  allow  for  better  comparisons.  Due  to  the  time  step 
being  derived  by  the  simulation  the  times  for  each  figure  could  not  perfectly 
match  but  the  physics  still  shows  the  same  situations  occurring  in  comparison 
to  Figure  38.  Image  (c)  for  all  figures  shows  the  initial  burst  as  the  two 
pressure  waves  meet  in  the  middle  of  the  simulation.  This  show  the  core 
radius  expanding  at  the  point  of  impact.  Image  (d)  for  all  the  simulations 
then  show  the  beginning  of  the  helical  instability.  This  shows  that  the  sensor 
allowed  for  the  natural  creation  of  the  instability  without  overpowering  it. 
As  the  number  of  cells  are  increased  in  the  z  direction  the  instability  becomes 
more  clear  as  can  be  seen  from  Figures  35  and  36  when  compared  to  Figures 
33,  34,  and  37.  The  increased  polynomial  order  did  not  show  the  helical 
instability  as  quick  as  the  other  examples  due  to  a  lower  resolution  in  the  z 
direction.  While  image (e)  in  the  figures  all  show  the  burst  cloud  becoming 
its  own  entity  in  the  center  of  the  simulations.  The  helical  instability  can 
easily  be  seen  in  the  (e)  image  for  all  the  figures. 

This  has  all  been  based  on  a  visual  inspection  therefore  a  more  quantifi¬ 
able  method  was  used  to  back  up  the  analysis  from  the  figures.  Since  this  is 
a  problem  based  on  the  Euler  Equations  and  the  boundary  conditions  have 
no  flow  passing  through  them  the  total  kinetic  energy  should  be  a  constant 
throughout  time.  Since  artificial  viscosity  is  being  applied  this  will  not  hold 
true  completely  and  therefore  if  too  much  viscosity  is  added  the  quality  of 
the  solution  decreases  and  it  can  be  tracked.  Figures  39  and  40  show  the 
total  kinetic  energy  of  the  problem  either  normalized  by  the  kinetic  energy 
of  the  initial  condition  or  the  error  when  compared  to  the  initial  condition. 
These  figures  both  show  that  there  is  only  small  incremental  improvement 
with  the  increase  of  cells  in  the  axial  direction  which  with  the  increase  in 
the  time  required  for  simulation  makes  a  solution  with  around  forty  cells  in 
the  axial  direction  seeming  to  be  most  optimal.  The  best  error  was  shown  to 
occur  on  the  increase  in  both  the  axial  and  radial  cells  which  lines  up  with 
the  other  charts. 

These  vorticity  magnitudes  and  kinetic  energy  plots  show  that  this  robust 
adaptation  of  the  sensor  was  able  to  approximate  the  same  viscous  physics 
effects  as  both  a  sixth  order  LES  and  DNS  simulations  without  requiring  the 
solving  of  the  RANS  equations  and  therefore  without  the  need  to  utilize  a 
turbulence  model.  These  figures  also  show  that  increasing  the  axial  resolution 
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has  only  a  slight  effect  in  the  decrease  of  the  error.  This  shows  the  superior 
results  to  be  for  the  9th  order  polynomial  with  a  grid  of  Nt  =  Nj  =  28,  Nk  = 
40  because  it  is  not  only  the  most  accurate. 
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(c)  t=4.0159 


(d)  t=4.5286 
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(e)  t=6.3687 


Figure  33:  Isosurfaces  of  vorticity  magnitude  for  a  polynomial  of  the  9th  order 
with  a  grid  Nt=Nj= 25  and  Nk= 30 
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(a)  t=0. 


(b)  t=1.8384 


(c)  t=3.8821  (d)  t=4.4949 
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(e)  t=6.3687 


Figure  34:  Isosurfaces  of  vorticity  magnitude  for  a  polynomial  of  the  9th  order 
with  a  grid  Nt=Nj= 25  and  Nk= 40 
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(a)  t=0. 


(b)  t=1.9658 


(e)  t=6.3687 


Figure  35:  Isosurfaces  of  vorticity  magnitude  for  a  polynomial  of  the  9th  order 
with  a  grid  Nt=Nj= 25  and  Nk= 40 
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(a)  t=0.  (b)  t=1.9064 


(e)  t=6.3687 


Figure  36:  Isosurfaces  of  vorticity  magnitude  for  a  polynomial  of  the  9th  order 
with  a  grid  Nt=Nj= 25  and  Nk= 45 
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(c)  t=3.9798 


(d)  t=4.5289 
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(e)  t=6.3687 

Figure  37:  Isosurfaces  of  vorticity  magnitude  for  a  polynomial  of  the  13th 
order  with  a  grid  Ni=N j=ll  and  N^= 25 
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(a)  t/Trot  =  1.01 


(b)  t/Trot  =  2.03 


(c)  t/Trot  =  3.05  (d)  t/Tmt  =  4.06 


Figure  38:  Isosurfaces  of  vorticity  magnitude  of  run  DNS2  (Trot  = 
^rc/V0max).[16\ 
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Figure  39:  Plot  of  the  integrated  total  kinetic  energy  error  versus  time  for 
the  cases 


Figure  40:  Plot  of  the  integrated  total  kinetic  energy  normalized  for  each 
case  versus  time 
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7  Conclusions 


In  conclusion  these  simulations  show  the  high  applicability  of  the  DG  method 
utilizing  the  Euler  equations  to  a  range  of  problems.  The  method  was  further 
improved  by  taking  the  Klockner  sensor  and  augmenting  it  to  improve  it’s 
capabilities  and  make  it  more  robust.  In  devising  a  new  system  with  which 
to  apply  the  baseline  modal  decay  to  the  sensor  allowing  for  the  sensor  to 
work  on  values  with  a  norm  value  of  zero  as  well  as  speed  which  was  one  of 
the  main  goals  of  this  report  and  its  most  important  findings.  The  savings 
in  cost  also  reduce  the  flop  count  in  a  two  dimensional  simulation  by  N'jp/ 2 
per  cell  per  timestep  and  for  three  dimensional  to  reduce  it  by  Nj  / 3  flops. 

This  sensor  was  then  taken  and  applied  to  the  three  dimensional  vortex 
bursting  problem  and  provided  accurate  results  in  showing  the  vortex  burst 
occurring  as  well  as  the  helical  instability  developing.  This  was  done  without 
the  need  to  work  with  either  a  LES  or  DNS  simulation  saving  immensely 
in  the  complexity  of  the  code  as  well  as  the  required  number  of  arbitrary 
coefficients  that  must  be  utilized  in  a  turbulence  model.  The  number  of  cells 
in  the  axial  direction  was  shown  to  be  the  most  important  factor  in  seeing  the 
development  of  the  helical  instability  while  the  resolution  in  the  transverse 
directions  played  more  of  a  role  in  the  pressure  distribution  charts.  The  9th 
order  polynomial  was  superior  to  the  13th  order  due  to  the  increase  in  speed 
it  allowed  as  well  as  the  better  overall  accuracy. 
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